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I. INTRODUCTION 



Water is one of the most studied substances in all its phases, vapor, liquid and solid - 
icei — due to its ubiquity in nature and great relevance to mankind. Despite the apparent 
simplicity of this molecule, it shows complex behavior and some of its properties remain 
poorly understood. Water's hydrogen-bonding and proton-disorder effects lead to a complex 
phase diagram, which has been progressively extended over many years. An extensive range 
of crystalline solid phases — or ice polymorphs — are known, most of which are stable and/or 
metastable under extreme conditions. On the other hand, different amorphous solid phases 
of water, or polyamorphs, have been discovered, including some ices which are identified as 
the most common water phases in the universe, being those found in interstellar space. 

Ice polyamorphs are usually distinguished by their characteristic densities. A low-density 
amorphous ice (LDA) was first synthesized in the 1930s by physical vapor deposition 2 and, 
more recently, by fast cooling of liquid water 3 . The existence of a second amorphous solid 
phase of high density (HDA) was established by Mishima and co-workers- in their exper- 
iments on the abrupt pressure-induced amorphization of hexagonal ice (Ih; the crystalline 
solid phase stable under ambient conditions) at low temperatures. After this discovery, 
many experimental and theoretical research efforts have been addressed to the characteriza- 
tion of the structural transitions between different crystalline and amorphous solid phases 
of water-. Diverse studies suggested the existence of two different mechanisms for the first- 
order pressure-induced transformation of Ih into HDA: at moderately low temperatures, the 
amorphization takes place by means of an endothermic melting of the crystalline morphol- 
ogy, whereas at very low temperatures the amorphous phase is the result of an exothermic 
structural collapse of the crystal^-. A reversible pressure-induced transformation of LDA 
into HDA was also first announced by Mishima^ and subsequently investigated in numer- 
ous experimental^"— and theoretical^"— studies. More recently, the existence of another 
amorphous solid phase with a very high density (VHDA) has been uncovered^ 1 ^^ and has 
become the subject of much research 2 ^"—. 

Studies on amorphous solid water have played a central role in the great interest on 
the understanding of the phenomenon of polyamorphism — the existence of different well- 
defined amorphous phases of the same substance — that has arisen in recent years 2 ^— . Solid 
water polyamorphism has stimulated the search for its liquid counterpart, i.e., the search 
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for two different liquids separated by a first-order phase transition and an associated liquid- 
liquid critical point^22r— . In this context, a extremely simple physical mechanism has been 
proposed as a generic explanation for the phenomenon of solid and liquid polyamorphism: 
the existence of a double well — or, more generically, of two characteristic length scales - 
in the intermolecular potential of a polyamorphic substance^™— . 

The extraordinary complexity of water's behavior has favored the development of a myr- 
iad of models intended for the study of distinct specific properties. Investigations on water 
polyamorphism in particular have taken great advantage of computer simulations based on 
multiple water models with very different levels of detail. For instance, there exist nu- 
merous studies on pressure- induced amorphization and other related amorphous transitions 
performed by molecular-dynamics simulations with atomistic water models' - —. However, 
an adequate understanding of the essential mechanisms of water polyamorphism — like the 
validity of the two-length-scales hypothesis — may require the exploration of more sim- 
ple or even minimal modeling approaches. Simple water models have been used for many 
years to study, for instance, its numerous thermodynamic anomalies and its structural order, 
phase diagram or solvation properties'^-. Some simple models for water polyamorphism 
have also been developed'^. On the other hand, there exist many simple generic models 
to test the validity of the two-length-scales hypothesis, mostly based on isotropic central 
potentials' 1 ^ - —. Models with anisotropic interactions are far more scarce due to the added 
complexity imposed by the directional bonds'"— . In the case of water, however, the strong 
impact of the directionality of the hydrogen bond on its properties makes the use of isotropic 
potentials particularly challenging, imposing a fine tuning of the model parameters in order 
to reproduce the desired properties" 

Perhaps the simplest model of water that incorporates a directional bonding scheme was 
introduced by Ben-Naim in the early 1970s to obtain a qualitative representation of the 
open hydrogen-bonded network of molecules that makes up liquid water— The model 
represents water molecules as two-dimensional Lennard- Jones disks with three equivalent 
hydrogen-bonding arms disposed at 120 degrees, as shown in Figured] The similarity of the 
shape of these simplified water molecules with a well known brand logotype has led to the 
adoption of the name "Mercedes-Benz" (MB) for the model. 

Despite its simplicity, Ben-Naim's MB model and its successive extensions and improve- 
ments have been shown to reproduce qualitatively different properties of water, includ- 
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ing some of its anomalies and the thermodynamic behavior of the melting transition—"—. 
Mercedes-Benz models have been used to study solvation and hydrophobicity problematic— 
and the properties of water under confinement^. Due to its flexibility, this class of models 
has been the subject of extensive analytical studies^"—. On the other hand and to the best 
of our knowledge, other more challenging characteristics of water, like the long disputed 
existence of two liquid phases or the properties of the solid amorphous phases and transi- 
tions, have not been explored to date in the minimal MB model. Regarding the main goal 
of our study, we consider that the MB model may be a useful anisotropic minimal modeling 
approach to study the essential mechanisms of water polyamorphism. In addition, bond- 
bending forces play a key role in many low- dimensional systems, such as in the amorphous 
freezing of soft polymer coils or silica nanoparticles in Langmuir monolayers^. In water, 
the interplay between the highly directional hydrogen-bonding network and the geometrical 
constraints determines the structure of liquid water and ice in two-dimensional layers, ei- 
ther under confinement or at open interfaces". Within this context, most computational 
studies have been devoted so far to the liquid structures", disregarding the behavior of 
amorphous solid phases. 

In summary, in this work we study for the first time the amorphous solid phases of the 
two-dimensional MB model of water and their transitions. In particular, we search for the 
existence of low-density (LDA) and high-density (HDA) solid amorphs and the determination 
of the nature of the transition between either hexagonal ice (Ih) and HDA as well as between 
LDA and HDA. Additionally, this approach allow us to study the validity of the two length 
scales hypothesis when directional bonds and a low dimensionality are introduced in the 
system. We place our results in the context of the known experimental results about ice. 

II. SIMULATION MODEL 

The MB pair-interaction potential is expressed as the sum of a radial and a directional 
term, 



U M b (n, fj) = U LJ ijij) + £/" H B {n, fj) . 



(1) 



The radial term, {7lj, is simply a Lennard- Jones potential, 




(2) 
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FIG. 1. Schematic representation of the two-dimensional Mercedes-Benz model of water (left) and 
its two crystalline morphologies or polymorphs (right). The water molecules are represented as 
Lennard- Jones disks with radius tlj combined with three hydrogen-bonding arms of length thb, 
here depicted as arrows. 

with the distance between the centers of the molecules as argument, = \fi — fj\. The 
directional term, {7hb; represents the water hydrogen bond and is denned by means of 
unnormalized gaussian functions, G (r) = exp [— r 2 /2a^ B ]. In his original model, Ben-Nairn 
defined Uub to be 

UBB(n, fj) = e H B G (nj - thb) B (fi, fj) , (3) 

where enB and thb are the depth and position of the bonding potential minimum, respec- 
tively, and 
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B (fi, fj) =J2 G (** ' ~ l ) G Oi ■ «w + 1) • ( 4 ) 

k,l=l 

Here it and ji are unitary vectors in the direction of every hydrogen-bonding arm of molecules 
i and j, respectively, whereas itij = fij/\fij\ is the unitary displacement vector between their 
centers. 

More recently, Silverstein and co-workers^ proposed a computational simplification of 
the model, by replacing expression fll]) by: 

B (fi, fj) = G (v(i, Uij) -1)G (w(j, itij) + 1) , (5) 

where 

v(i, Uij) = max(?i • , i 2 ■ u^, z 3 • Uij) (6) 

w(j, u^) = min(j! • uy, % ■ 0,^,% ■ u^) (7) 

According to the previous definitions, the MB model has two different bonding distances 
given by the minimum of the Lennard- Jones potential, r^j = 2 1//6 <7lj, and the hydrogen bond 
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length, thb, introduced in Eq. d3J). As a consequence of these two bonding mechanisms, two 
crystalline solid morphologies can be found in the model: a low-density hydrogen-bonded 
honeycomb lattice and a high-density triangular lattice of Lennard- Jones disks, as shown in 
Figure [TJ In particular, it has been shown by means of Monte Carlo NPT simulations that 
the melting of the MB honeycomb structure reproduces qualitatively the thermodynamic 
properties of the melting of ice Ih^i. Therefore, it is reasonable to expect the existence in 
the model of at least two solid amorphs with low and high characteristic densities, which 
eventually could be associated with the low- and high-density amorphous ices. 

In order to explore the existence and characteristics of solid amorphs in the two- 
dimensional MB model, we performed extensive equilibrium Monte Carlo NPT simulations 
with a system composed with up to 1200 MB particles in a rectangular cell with periodic 
boundary conditions. At least 50 independent runs of 2 ■ 10 7 steps were performed for every 
point using the model parameters proposed in previous works^: €lj = 0.1, cxlj = 0.7, 
£hb = 1-0, &hb = 0.085, r HB = 1.0. After equilibration, measures of the internal energy, 
volume and structure were taken and averaged over all runs. As usual, the system enthalpy 
and heat capacity were calculated as: 



Here (. . .) denotes averages over runs and the parameters are expressed in reduced units, 
relative to the hydrogen bond parameters: T* = k B T j\e HB \, V* = V/r 2 HB) U* = U/\€hb\, 



In the next section we present the results obtained from our simulations. In many cases, 
they correspond to simulations performed under very low temperature conditions. Under 
such circumstances one must be aware of the difficulties of obtaining well equilibrated struc- 
tures at the transition region using simple NPT Monte Carlo simulations; thus extensive 
computer work is required to avoid the system being trapped into a local minimum. In ad- 
dition, the simplicity of the MB model — as with any minimal model — makes comparison 
with the experiments relevant only on a qualitative level. 




(8) 
(9) 



H* = H/\e HB \,P* = r 2 HB P/\e HB 
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III. RESULTS AND DISCUSSION 



In our simulations, we tried to reproduce the experimental amor phizat ion paths estab- 
lished by Mishima and collaborators in their pioneering works on amorphous ices. In partic- 
ular, we focus on the amorphization of Ih ice — which we identify with the honeycomb lattice 
- into HDA by compression at very low temperature^ and on the reversible transformation 
between LDA and HDA ices obtained by compression and decompression with annealing^. 
Except for the latter case, we worked well below the melting point of the honeycomb crystal, 
~ 0.15 at P* = 0.1 61 , assuming that the resulting sample structures remain in a solid 
state or, at least, in a very viscous amorphous phase. We shall discuss this assumption on 
the basis of the rigidity percolation theory applied to amorphous solids. 



A. Pressure-induced amorphization of ice Ih 

The transformation of ice Ih into HDA is studied by compressing a sample of N molecules, 
disposed initially in the honeycomb lattice, at a very low temperature, T*=0.05. Figure 
[2] provides a first insight into the general behavior of the system during this process for 
N = 1200 MB particles. As the pressure is increased, the system responds initially with 
just a slight reduction of the volume and a small displacement of the particles from their 
equilibrium positions, while keeping the global honeycomb structure. The typical crystalline 



morphology at P* = 0.6 is shown in the inset of Fig. 2(a) Consistently, the radial distri- 



bution function at P* = 0.6 (see Fig. 2(b)) identifies two clear maxima corresponding to 
the first and second nearest neighbor positions in the compressed honeycomb lattice. At 
around P* = 0.7, an abrupt collapse of the stressed honeycomb structure takes place, lead- 
ing to the arrangement of the particles into a high-density amorphous configuration, which 
we identify with HDA ice. As revealed by its radial distribution function at P* = 1.0 (see 



Fig. 2(b) ), this amorphous structure is associated with a significant formation of LJ bonds, 
corresponding to the peak at around the LJ equilibrium distance, r^j 0.78, which replace 
a fraction of the original HB bonds of the honeycomb lattice, indicated by the peaks at 
around thb = 1.0. A snapshot of HDA ice at P* = 1.0 is also shown in the inset plot of Fig. 



2(a) to compare with the crystalline structure. The HDA morphology remains with little 



change when the system is further subjected to an isothermal decompression. Qualitatively, 
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FIG. 2. (a) Pressure-induced amorphization of an Ih crystal sample with N = 1200 at T*=0.05 
(upper curve) and decompression of the resulting amorph at the same temperature (lower curve), 
with insets showing the morphologies found at low and high pressures for the compression curve, 
(b) From the same sample, radial distribution functions for the stressed Ih crystal (P* = 0.6) and 
for the high-density amorphous phase (P* = 1.0). (c) Probability histograms of the conflgurational 
energy for three selected pressures of the compression process. 



this behavior of the system volume is completely consistent with what can be observed in 
experiments^ and is a clear indication of a pressure-induced phase transition, probably of a 
first-order kind, as shown by the abrupt drop in the volume even for such a relatively small 
system. 

We have further investigated the nature of the Ih— )-HDA transition by studying different 
parameters. The analysis of the probability distribution function of the total conflgurational 
energy (Fig. 2(c)) clearly identifies single peaks before (P* = 0.68) and after (P* = 0.78) 



the transition, whereas for pressures close to the transition point (P* = 0.73) two maxima 
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FIG. 3. Upper panel: evolution of the mean number of HB and LJ bonds per particle, n&, along the 
Ih— 7-HDA transition for N = 1200; see the text for the bond-counting criterion used. Lower panel: 
corresponding evolution of the system enthalpy, H*/N, with the projection and intersection of the 
lines from each side of the transition used to estimate the thermodynamic transition pressure, Pq. 

are found. This result is a strong indication of the first-order nature of the transitional^. 

Another way to characterize this transition is by studying the evolution of the bonds 
within samples as the pressure is increased. Qualitatively, it is evident that the structure 
must evolve from a rigid honeycomb lattice, connected by just HB bonds, into a completely 
different rigid structure, presumably independent from the former in the limit of very high 
pressures, composed of a triangular lattice of particles in close contact and governed by the 
soft-core barriers of the LJ potential. In order to obtain some further insight into how this 
evolution takes place, we computed for every pressure the mean number of bonds of every 
type and the connectivity of the network denned by all the bonds. The criterion used to take 
bonds into account has been the following: a bond, either of LJ or HB type, is considered 
as established between any two given particles when the strength of the interaction is above 
0.75 of its maximum possible value. In the case of the LJ potential, since it represents a 
soft-core barrier for high pressure configurations, we also consider the bond established when 
the interparticle distance is below its optimum value, tlj. The upper panel of Figure [3] shows 
the evolution for the split mean number of HB and LJ bonds obtained for A" = 1200. As 
expected, there is a sigmoidal-shaped increment of the number of LJ bonds and a reduction 
of the number of HB bonds in the transition region. The total number of bonds of any 
kind increases monotonically as one would expect from the different maximum coordination 
number of both lattices. However, it is remarkable that the number of HB bonds decreases 
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very slowly after the transition region and remains significant at relatively high pressures, 
indicating that some HB bonds still exist within compact configurations. This behavior 
has an impact on the thermodynamic properties of the transition. The lower panel of 
Figure 13] shows the system enthalpy per particle for N = 1200 as a function of the pressure. 
The notable step down shown by the enthalpy at the transition region is numerically a 
consequence of the significant persistence of HB bonds after the collapse of the honeycomb 
structure, leading to a relatively small increase in the internal energy, AU*/N w 0.1, in 
front of the considerable decrease of the system volume, AV*/N rs 0.5. Thermodynamically, 
such drop of enthalpy is a clear indication of the release of a hysteresis heat corresponding 
to the system relaxation from a metastable state: since the temperature is so low, the 
system gets kinetically trapped into the crystal phase until the overpressurization is high 
enough to overcome the energy barriers. Therefore, the thermodynamic transition point 
can be estimated from the intersection of the projected enthalpy lines from each side of the 
transition, as shown in Figure EJ From this calculation we get a value for the transition 
pressure of P * = 0.43 ± 0.07. 

Finally, the identification of the bonds allows us to study the clustering of the networks 
of bonded particles. In all cases we found that the connectivity of the network of bonds is 
maintained during the transition, so that all the particles remain connected into a single 
cluster at all pressures. According to the rigidity percolation theory^ 9 -, this behavior - 
in combination with the monotonic increase in the mean number of bonds — indicates 
that the solidity of the sample structure is maintained during the transition. Therefore, 
this suggests that amorphization occurs via mechanical collapse instead of a melting of the 
crystal structure. 

All these observations are consistent with the known experimental and simulation results 
on the pressure-induced amorphization of ice Ih at very low temperatures^ - -. 

B. Transformations between LDA and HDA ices 

The second process explored in our simulations with the MB model is the reversible 
transformation between high- and low-density amorphous structures. As in the previous 
case, we apply a procedure equivalent to Mishima's experiments to simulate a low-density 
amorphous solid, or LDA ice, and its reversible transformation into HDA 9 . Specifically, the 
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transformation HDA— s-LDA has been obtained by applying decompression with annealing, 
i.e., by increasing the temperature of the HDA sample as the pressure is lowered, whereas 
the reverse transformation HDA— t-LDA has been achieved by compressing the LDA ice at 
high pressure under very low temperature conditions. 




(a) 





(b) (c) 

FIG. 4. Results for the reversible transformation between low-density (LDA) and high-density 
(HDA) amorphous structures by compression and decompression with annealing, (a) Evolution 
of the system specific volume for the compression at T* = 0.06 (upper curve of the lower panel) 
and decompression with annealing (lower curve of the same panel). The annealing consists of a 
linear increase of the temperature from T* = 0.06 to T* = 0.13 (upper panel), (b) Detail of the 
LDA morphology, (c) Radial distribution functions of the corresponding high-density (HDA) and 
low-density (LDA) amorphous phases. 



Figure H] shows the main results obtained from the indicated transformation procedures 
on a sample of size N = 1200. First, a closed transformation cycle LDAt>HDA has been 



successfully achieved, as shown in the lower panel of Figure 4(a) The LDA structure has 
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been produced by a linear increase of the system temperature from 0.05 to 0.13 as the pres- 



sure was reduced from 1.2 to 0.01 (upper panel of Figure 4(a) ). The mean number density of 
the resulting structure, which is mainly controlled by the final temperature, is Plda ~ 0.71, 
a value slightly lower than that corresponding to the honeycomb lattice, p\ h ~ 0.77. Its 



radial distribution function, shown in Figure 4(c) , confirms that the LDA morphology is the 



amorphous counterpart of the honeycomb lattice, being mainly composed of HB bonds but 
with many structural defects. This morphology remains intact when the temperature is set 
back to T* = 0.06. By applying an increasing pressure at such a low temperature, the LDA 
morphology experiences a gradual compaction to arrive once more at the HDA structure. 
The continuous, smooth nature of the transformation between these amorphous forms is 
not what is found in experiments and simulations with other more realistic water models, 
from which it has been well established that its true nature is that of a first-order phase 
transition, with associated latent heats*^. We tested also other configurations of poten- 
tial LDA structures with a somewhat higher density, produced by reducing the maximum 
temperature of the annealing process. In all cases — including some with a density even 
slightly higher than that corresponding to the honeycomb lattice — the same qualitative 
results were obtained. For higher maximum annealing temperatures, a complete melting of 
the structure is soon obtained. Therefore, we were unable to find in the two-dimensional 
MB model any low-density amorphous structure with enough structural stability to produce 
a pressure-induced discontinuous phase transition into a high-density amorphous form. In- 



deed, we want to stress that the transformation cycle shown in Figure 4(a) corresponds just 
to the results qualitatively closer to Mishima's experiments that we were able to find in our 
simulations. In particular, the shape of the decompression HDA— ^LDA curve is controlled 
by the annealing temperatures and therefore can be strongly distorted by using different 
annealing conditions. 

We consider that the origin of the apparent mechanical instability of the low-density 
amorphous phase in this model is most probably related to the low maximum coordination 
number imposed by the HB bonds and its interplay with the low dimensionality of the 
system, which geometrically forbids the existence of defects, inherent to any amorphous 
structure, without an associated reduction of the mean number of directed bonds. As can 



be observed in the example of Figure 4(b) , most defects of the LDA structure are associated 



with misalignments of the directed bonding arms. These misalignments have effects at scales 
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larger than the distance of first-nearest neighbors: as can be observed, the formation of non 
hexagonal cells — closed loops of either less or more than six elements — is very frequent. 
This effect, combined with the limited possibilities of tessellation of the two-dimensional 
space, imposes the existence of many unbonded arms. For instance, in the case illustrated 
by Figure HJ the total mean coordination number - - calculated by means of the bond- 
counting criterion introduced in the previous section — is 2.70, or just 2.43 if only the HB 
bonds are taken into account. Obviously, any significant decrease of the mean number of HB 
bonds in this model implies a considerable increase in the total configurational energy of the 
sample: continuing with the example from Figure IU the mean configurational energy per 
particle of such an LDA structure is approximately -1.34, almost 15% higher than the energy 
corresponding to the unstressed honeycomb lattice, -1.57. Such an energy is still significantly 
higher than that of the stressed honeycomb lattice at the point of collapse, -1.45, showing 
the overall weakness of the structure. This point represents a significant difference with 
respect to what is observed in three-dimensional simulations with tetrahedral water models, 
in which LDA ice keeps the fully coordinated network structure^. 

IV. CONCLUDING REMARKS 

We have performed extensive NPT Monte Carlo simulations of the two-dimensional MB 
model in order to study the essential underlying physical mechanisms behind ice polyamor- 
phism. In particular, we have investigated the validity of the two-length-scales hypothesis, 
previously suggested as the minimal ingredient for the interaction potential of polyamorphic 
materials, when directional bonds and a low system dimensionality are considered. 

To this end we have investigated, in the first place, the pressure-induced transforma- 
tion of ice Ih into HDA at very low temperatures. Our results suggest the existence of a 
first-order phase transition in which amorphization occurs via mechanical collapse of the 
crystal honeycomb lattice from a kinetically trapped metastable state into HDA ice. This 
amorphous structure is associated with a significant formation of LJ bonds that replace a 
small fraction of HB bonds in the original crystal. This mechanism ensures the network 
connectivity during the transition, thus preventing the system from melting. This result is 
in agreement with the experimental observations of pressure-induced amorphization of ice 
Ih under very low temperature conditions 4,7 . 
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In the second place we have explored the transformation between high- and low-density 
amorphous ices by performing an (isothermal compression)-(annealed decompression) cycle. 
Our results indicate the existence of a continuous transformation between such amorphous 
structures that is in contradiction with the experimental findings 9,80 . We consider that this 
discrepancy can be attributed to the low coordination of the low-density amorphous phase, 
which has no significant structural stability. This low connectivity arises as a consequence of 
the constraints imposed by the bond directionality and the low dimensionality of the system. 
Therefore our results provide a clear indication that an effective interaction potential with 
two characteristic length scales does not guarantee by itself a first-order phase transition 
between polyamorphs when it is accompanied by strong bonding constraints. 

We hope that these results might stimulate new experiments performed in low dimensional 
systems to study the effect of geometrical constraints and the validity of the predictions of 
the minimal MB model. 
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